#! /usr/bin/python
# -*- coding: utf-8 -*-

# create_coeff_sequences.py
import math
import random

#approximation of the expected value of 1/x, where x is a random
#variable distributed evenly between 0 and 1.
#To compute this approximation, I subdivide [0,1] into n subintervals
#and compute 1/p for the midpoint p of each subinterval. These are
#summed up and divided by n. 
def one_by_x(n):
    sum = 0
    point = 1 / (2 * n)
    width = 1 / n
    for i in range(n):
        sum = sum + 1 / point
        point = point + width
    return sum / n

#It turns out that the above function is roughly the natural logorithm shifted
#by 2. I guess this can easily be explained considering that 1/x is the
#derivative of ln x. 
def one_by_x_via_log(n):
    return 2 + math.log(n)

#rand_factor computes a random factor between 1 and inf. The bigger you
#choose factorW, the bigger the factor. 0.1 is a value for factorW
#that seems to give nice results. 
#def rand_factor(factorW):
#    return 1 + factorW / random.random()
def rand_factor(factorW):
    return pow(random.random(), -factorW)

#Based on the same calls of rand_factor, I compute 2 sequences: one
#where I create "float coefficients" that are floored "at the last
#moment" (= smooth), and one where the new coefficient is immediately
#floored after multiplication with factor. The first approach gives a
#smoother growth of coefficients whereas the second generates more
#identical coefficients because it is less likely to increase the
#current coefficient.   
def both_coeff_seq(m,floatCoeff,factorW):
    smoothCoeff = floatCoeff
    coeff = math.floor(floatCoeff)
    smooth_sym = 0
    sym = 0
    for i in range(m-1):
        print(coeff, math.floor(smoothCoeff))
        factor = rand_factor(factorW)
        if math.floor(smoothCoeff) == math.floor(smoothCoeff * factor):
            smooth_sym = smooth_sym + 1
        smoothCoeff = smoothCoeff * factor
        if math.floor(coeff * factor) == coeff:
            sym = sym + 1
        coeff = math.floor(coeff * factor)
    print(coeff, math.floor(smoothCoeff))
    print("Symmetries:")
    print(sym,smooth_sym)

#Now the smooth sequence only: 
# def smooth_coeff_seq(m,floatCoeff,factorW):
#     smoothCoeff = floatCoeff
#     smooth_sym = 0
#     oldIntCoeff = 0
#     for i in range(m):
#         factor = rand_factor(factorW)
#         smoothCoeff = smoothCoeff * factor
#         newIntCoeff = math.floor(smoothCoeff)
#         if i==0:
#             smallest = newIntCoeff
#         else:
#             if newIntCoeff == oldIntCoeff:
#                 smooth_sym = smooth_sym + 1
#         oldIntCoeff = newIntCoeff
#         print(newIntCoeff, end="\;")
#     if (m>1):
#         print("Excentricity: ", pow(newIntCoeff/smallest, 1/(m-1)))
#         print("Symmetry quotient: ",smooth_sym/(m-1))

def smooth_coeff_seq_latex(m,floatCoeff,factorW):
    print(factorW, end=" & ")
    smoothCoeff = floatCoeff
    smooth_sym = 0
    oldIntCoeff = 0
    intCoeff = [0] * m
    for i in range(m):
        factor = rand_factor(factorW)
        smoothCoeff = smoothCoeff * factor
        intCoeff[(m-1)-i] = math.floor(smoothCoeff)
        if i>0 and intCoeff[(m-1)-i] == intCoeff[m-i]:
             smooth_sym = smooth_sym + 1
    for i in range(m):
        print(intCoeff[i], end="\;")
    if (m>1):
        print(" & ", pow(intCoeff[0]/intCoeff[m-1], 1/(m-1)), end="")
        print(" & ",smooth_sym/(m-1),end="")
    t1 = intCoeff[0]
    t2 = 0
    for i in range(math.floor(m/2),m):
        t2 = t2 + intCoeff[i]
    t4 = 0
    for i in range(0,math.floor(m/2)):
        t4 = t4 +  intCoeff[i]
    t3 = math.ceil((t2 + t4) / 2)
    t5 = t2 + t4 - intCoeff[0]
    print(" & ",t1,"\;",t2,"\;",t3,"\;",t4,"\;",t5,"\\\\")
    
def smooth_coeff_seq(m,floatCoeff,factorW):
    smoothCoeff = floatCoeff
    oldIntCoeff = 0
    intCoeff = [0] * m
    t = [0] * 5
    for i in range(m):
        factor = rand_factor(factorW)
        smoothCoeff = smoothCoeff * factor
        intCoeff[(m-1)-i] = math.floor(smoothCoeff)
    t[0] = intCoeff[0]
    t[1] = 0
    for i in range(math.floor(m/2),m):
        t[1] = t[1] + intCoeff[i]
    t[3] = 0
    for i in range(0,math.floor(m/2)):
        t[3] = t[3] +  intCoeff[i]
    t[2] = math.ceil((t[1] + t[3]) / 2)
    t[4] = t[1] + t[3] - intCoeff[0]
    for j in range(5): 
        print(m,end=" "),
        print(t[j],end=" ")
        for i in range(m):
            print(intCoeff[i],end=" ")
        print("")

def smooth_coeff_seq_collection(m,floatCoeff,factorW, I):
    for i in range(I):
        smooth_coeff_seq(m,floatCoeff,factorW)

def iwoca_collection(m,I):        
    smooth_coeff_seq_collection(m,1.0,0.1, I)
    smooth_coeff_seq_collection(m,1.0,0.2, I)
    smooth_coeff_seq_collection(m,1.0,0.3, I)
    smooth_coeff_seq_collection(m,2.0,0.1, I)
    smooth_coeff_seq_collection(m,2.0,0.2, I)
    smooth_coeff_seq_collection(m,2.0,0.3, I)
  
